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FORTRAN IV COMPUTER PROGRAM FOR CALCULATING CRITICAL 
SPEEDS OF ROTATING SHAFTS 
by Roger J. Trivisonno 
Lewis Research Center 

SUMMARY 

This report describes a FORTRAN IV computer program, written for the IBM DCS 
7094/7044 computer, that calculates the critical speeds of rotating shafts. The shaft 
may include bearings, couplings, extra masses (nonshaft mass), and disks for the gyro- 
scopic effect. Shear deflection is also taken into account, and provision is made in the 
program for sections of the shaft that are tapered. The boundary conditions at the ends 
of the shaft can be fixed (deflection and slope equal to zero) or free (shear and moment 
equal to zero). The fixed end condition enables the program to calculate the natural fre- 
quencies of cantilever beams. Instead of using the lumped -parameter method, the pro- 
gram uses continuous integration of the differential equations of beam flexure across dif- 
ferent shaft sections. The advantages of this method over the usual lumped -parameter 
method are less data preparation and better approximation of the distribution of the mass 
of the shaft. 

A main feature of the program is the nature of the output. The Calcomp plotter is 
used to produce a drawing of the shaft with superimposed deflection curves at the critical 
speeds, together with all pertinent information related to the shaft. 

INTRODUCTION 

Most methods of calculating critical speeds are based on the lumped -parameter 
method, such as that of Prohl (ref. 1). This lumped -parameter technique consists in 
transforming a shaft system into a series of mass points so that their spacing approx- 
imates the distribution of the mass in the actual shaft. The next step is to integrate from 
mass point to mass point for the purpose of determining the shear, moment, slope, and 
deflection at each mass point. Consequently, in order to represent accurately the mass 
distribution of a shaft system, a considerable amount of data preparation is necessary. 



This is especially true for sections of the shaft that are tapered, where it is necessary to 
approximate the taper with a series of small increments. 

The method described in this report improves on the lumped-parameter method by 
integrating along each shaft section by means of the fourth -order numerical integration 
technique of Runge-Kutta. The advantages of this method are that both the time for data 
preparation and the chance for input error are considerably reduced and that the distri- 
bution of the mass of the shaft is better approximated. 

The three main sections of this report are STATEMENT OF PROBLEM, COMPUTER 
PROGRAM, and the appendixes. The first section contains the statement of the problem 
along with the method of solution. The second section contains the program description, 
the FORTRAN symbols for the NAMELIST card, and the method in which a shaft system is 
broken down in sections and stations for proper input to the program. Also, this section 
presents example problems, together with printed and plotted output, demonstrating the 
flexibility of the program. Finally, the appendixes contain the definitions of symbols and 
the listing of the FORTRAN program. 


STATEMENT OF PROBLEM 


A shaft of variable diameter, mass density, and modulus of elasticity having bear- 
ings, couplings, disks, and extra masses (nonshaft mass) is defined to be a complex 
shaft system. When such a system is rotating steadily at a critical speed co, the follow- 
ing boundary value problem has a solution: 


_1_/ei5!iW 2 y 

dx 2 \ dx 2 / 


where 


2Y = ^ 2 y 

dx 


dM 

dx 


= V M = El 


d 2 Y 

dx 2 


dY 

dx 


= e 


2 

E(x) modulus of elasticity at x, lb/in. 

I(x) diametral moment of inertia at x, in.^ 

2 p 

?7(x) mass per unit length of a general cross section, lb -sec /in. 
Y(x) deflection at x 
M(x) moment at x 


(1) 
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V(x) shear to right of x 

0(x) slope at x 

x axial position on shaft 

All symbols are defined in appendix A. 

The terminal boundary conditions are introduced by (1) free overhang 

M(0) = V(0) - 0 


or (2) fixed end 


# 


Y(0) - 0(0) = 0 


and free overhang 

M(L) = V(L) = 0 

The interior boundary conditions are introduced by (1) rigid bearing at station n 

Y(x) = 0 
0,M continuous 

or (2) coupling at station n 

M(x) = 0 
Y,V continuous 

The interior discontinuities are introduced by (1) flexible bearing at station n 

V(x + ) = V(x“) - KY(x) 

(2) extra mass at station n 

V(x + ) = V(x') + S n co 2 Y(x) 

or (3) gyroscopic disk at station n 

M(x + ) = M(x') + (A - B)co 2 e(x) 
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where 


K bearing stiffness factor (spring rate), lb/in. 

o 

A mass moment of inertia of disk about its axis of symmetry, lb -in. -sec 

B mass moment of inertia of disk about axis through center of gravity and normal to 

2 

axis of symmetry, lb -in. -sec 

2 

S n nonshaft mass at station n, lb -sec /in. 

Method of Solution 

To obtain the critical speeds of a shaft system, the shaft is transformed into a se- 
ries of sections and stations. Stations are axial positions on the shaft where disks, 
bearings, couplings, and extra masses are located or where the formulas for E, rj, 
and I change. A section is the portion of the shaft between two stations. By using the 
following method on a shaft system broken down in this manner, the critical speeds can 
be found. 

The method can best be illustrated by studying a simple uniform shaft, with no inte- 
rior boundary conditions or discontinuities. The boundary conditions at both ends of the 
shaft are free. 

Equation (1) is composed of the following first-order linear differential equations of 
beam flexure: 


v(x) = 7)(x)cl)^Y(x) 

(2) 

— M(x) = V(x) 

(3) 

dx 

d e(x) = 
dx EI(x) 

(4) 

— y(x) = e(x) 

(5) 

dx 


To obtain a simultaneous solution for this system of equations, the fourth -order 
Runge-Kutta method is used. The boundary conditions at x = 0 are used for the initial 
conditions. This system of equations is linear and homogeneous. The theorem of the 
superposition principle states that any linear combination of two solutions of a linear 
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homogeneous differential equation is again a solution. Since Y(0) and 0(0) are unspec- 
ified, we break the problem down into the following two cases: 


Case 1 

Case 2 

o 

n 

i— i 

o 

> 

< 

o 

CO 

II 

o 

M 0 i = ° 

M 02 ~ 0 

e oi : 0 

°02 = 1 

Y 01 =1 

Y 02=° 


where the subscripts denote the station and case number, respectively. 

After the integration across the entire shaft is completed, the solution to the initial - 
value problem with Y(0 ) = Yq and 0(0 ) = 0q is 


V n - V nl Y 0 + V n2®0 

(6) 

M n = M nl Y 0 + V n2 e 0 

(7) 

e n = e nl Y 0 e + e n2 e 0 

(8) 

Y n = Y nl Y 0 + Y n2°0 

(9) 


Now by applying the boundary condition at x = L, V(L) = 0, in equation (6) and solving for 
0 q and substituting in equation (7), we have: 


M 


L3 




M n2 V nl\ 
V n2 / 


( 10 ) 


When an a) is used that makes equal zero, we have satisfied the four bound- 

ary conditions and the problem is solved. Hence, by plotting M-^, in equation (10), as 
a function of the assumed speed co and noting where becomes zero, the critical 

speeds are obtained. 
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Shear Deflection 


Equation (5), which defines deflection, considers bending stresses only. There is, 
however, some additional deflection due to shear. Usually this is negligible; however, 
in cylindrical shafts with small span-to-depth ratios, stresses are likely to be high. 
Consequently, the deflection due to shear constitutes a considerable part of the total de- 
flection and is taken into account by Roark (ref. 2) as follows: 

— Y (x) = F < x)V ( x) (11) 

dx s A(x)G(x) 

where 

G(x) = E(x) - 

2(1 + v) 

o % 

F(x) = (7 + M(l + Q + (20 + 12fc*)c 

2 ^ 

6(1 + i/)(i + o 

A(x) is the cross-sectional area at x, v is Poisson T s ratio (= 1/3), and c = D^/D q . 

Now subtracting equation (11) from equation (5) yields 

— Y(x) = 0(x) - y(x)V(x) (14) 

dx 


( 12 ) 

(ref. 3) (13) 


where 


y(x) = - 

A(x)G(x) 

Interior Boundary Conditions 

In the previous section a simple uniform shaft was used to illustrate the method of 
calculating critical speeds. When internal boundary conditions are added (i.e. , rigid 
bearings and couplings), creating a multispan shaft system, the method of calculating 
critical speeds changes only at the boundary conditions. The explanation of these 
changes follows. 
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Rigid bearing. - This type of bearing is infinitely stiff against displacement but of- 
fers no restraint to tilting of the shaft. At the point on the shaft where the rigid bearing 
is located, x = x R , the new boundary conditions are zero deflection, continuous 6 and 
M, and an unknown jump in V. 

Equations (6) to (9) define V, M, Q, and Y in the span from x = 0 to x = x R . 

Now when we apply the condition Y(x^) = 0 to equation (9), it follows that 

Y(0) = C R 0(O) (15) 


where 


p Y n2^ x R ) 

^ Y (x ) 

nl' x R 

and equations (6) to (9) simplify to 

V n3 = (V n2 + C R V nl )e 0 
M n3 = (M n2 + C R M nl )e 0 
e n3 = ^ e n2 + C R e nl )e 0 
Y n3 = (Y n2 + C R Y nl' e 0 


(16) 

(17) 

(18) 
(19) 


Since V(x R ) = V R is unspecified, the problem is broken down for the next span into 
the following two cases: 


Case 1 

Case 2 

V R1 =1 

V R2 = ° 

M R1 = ° 

M +R2 = M -R2 - X R^ + C R M 

e Rl = 0 

®+R2 = ®-R2^ x r' + C R e - 

Y R1=° 

Y R2 = 0 


7 



For points in this span the solution to the boundary value problem with boundary con- 
ditions e(0) = 6 q and V(x R ) = V R is given by 


= V nl V R + V n2 e 0 

(20) 

= M nl V R + M n2 e 0 

(21) 

= e nl V R + 0 n2 e 0 

(22) 

= Y nl V R + Y n2«0 

(23) 


Couplings . - At the point on the shaft where a coupling is located, x = x c , the new 
boundary conditions are zero moment, continuous Y and V, and an unknown jump in 0. 
Now when we apply the condition M(x^) = 0 to equation (7), it follows that 

Y(0) = C C 0(O) (24) 


where 



M 


n2 


M 


nl 


( x c ) 


and equations (6) to (9) simplify to 


V n3 = (V n2 + C C V nl' e 0 

(25) 

M n3 = (M n2 + C C M nl )e 0 

(26) 

e n3 = ^ e n2 + C C e nl )e 0 

(27) 

Y n3 = ^ Y n2 + C C Y nl )e 0 

(28) 


Since 0 (xq) = 0q is unspecified, the problem is broken down for the next span into the 
following two cases: 
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Case 1 


Case 2 


V C1 = 0 V +C2 = V_ c2 (x c ) + C c (x c ) 

M C1 = 0 M C2 = 0 

e cl = i e +c2 = 0 

Y c 1 = 0 Y +C2 = Y_ c2 (x c ) + C C Y c 1 (x c ) 

This completes the treatment on interior boundary conditions. 


Interior Discontinuities 

This section treats interior discontinuities introduced by flexible bearings, extra 
masses (nonshaft mass), and disks attached to the shaft. The flexible bearings and extra 
masses introduce a discontinuity in shear. A disk attached at some point on the shaft in- 
troduces to the rotating system a gyroscopic effect which causes a discontinuity in the 

moment. The treatment of these three cases follows. 

Flexible bearing. - This type of bearing is elastic against displacement but offers no 
restraint to tilting of the shaft. At thfe point on the shaft where the flexible bearing is lo- 
cated the shear is changed by subtracting from it the product of the spring rate K of 

’ ! . ! 
the bearing and the deflection Y(x) 

V(x + ) = V(x~) - KY(x) (29) 

Extra masses (nonshaft mass). - Masses that are attached to the shaft but do not 
contribute to the stiffness of the shaft bring about a change in shear at the point on the 
shaft where they are concentrated. By using the following equation, the shear to the 
right of the station where the mass is located can be computed: 

V(x + ) = So) 2 Y(x) + V(x”) (30) 

where S is the extra mass at x. 

Gyroscopic effect of disks attached to shaft. - The gyroscopic effect of disks at- 
tached to the shaft can be taken into account for calculating the critical speeds by using 
the following formula from reference 1: 
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M(x + ) = (A - B )w 2 e(x) + M(x”) 


(31) 


where 

A mass moment of inertia of disk about its axis of symmetry 

B mass moment of inertia of disk about axis through center of gravity and normal to 
axis of symmetry 

For a flat disk, A and B are given by 


7rphD d 

32 


rrphD^j 

B = 2 . 

64 




Deflection Curve Calculation 

An important phase of analyzing complex shaft systems is that of studying the deflec- 
tion curve at a critical speed. One of the important aspects of a deflection curve is that 
it reveals where the greatest bending of the shaft is taking place at a critical speed. 
Knowledge of this nature enables the engineer to make the necessary changes in the shaft 
design to meet the requirements of the design running speed. As an aid to the design 
engineer, the program superimposes all deflection curves at the critical speeds on a 
drawing of the shaft. 

By using the following formula, a value for Y n3 can be obtained for every x n : 

Y n3 = ( Y n2 + C R,C, F Y nl )e 0 ( 34) 

where is used from a boundary condition to a rigid bearing and is equal to 
"t Y n2^ x R )//Y nl^ x R )], is used from a boundary condition to a coupling and is equal to 

“f M n2^ x C^ M nl^ x C^’ and C F is used from a boundary condition to a free overhang and 
is equal to -[V n2 (x L )/V nl (x L )] . 

It should be noted that the deflection curve produced by equation (34) represents the 
shape of the deflected shaft centerline at the critical speed. It does not represent the 
true magnitude of the shaft deflection. 
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COMPUTER PROGRAM 


The critical -speed program is written in FORTRAN IV language for the IBM DCS 
7094/7044 computer. It is made up of a main program, six subroutines, and one library 
subroutine used for plotting, which is not included in this report. The program can solve 
a multispan shaft system for a maximum of 200 stations and any combination of bearings 
and couplings up to nine. See appendix B for the listing of the FORTRAN program. 

The program will compute the first three critical speeds in the rpm interval speci- 
fied. If more than three critical speeds are desired, another computer run is made with 
the rpm interval set at a higher level. 

The execution time to compute a critical speed depends, of course, on the shaft 
length and the number of rpm increments that will be used. For example, in the section 
on sample computer problems, problem 1 involves a shaft 53.39 inches long with 51 
shaft sections. For 150 increments of 200 rpm each, the execution time is approx- 
imately 1 minute. 

A more detailed explanation of the program follows. 


Program Description 

Main program CSP . - The main program performs the following functions: 

(1) Reads all input data 

(2) Writes all input data 

(3) Sets switches to regulate the flow of the program 

(4) Calculates gyroscopic data (if any) 

(5) Calls subroutine MAIN 

Subroutine MAIN . - This subroutine (1) sets up initial conditions for numerical in- 
tegration and then calls subroutine RUNGEK; (2) performs calculations for changes 
brought about by extra masses, flexible and rigid bearings, couplings, and gyroscopic 
effects; (3) calculates critical speed and writes converged value; (4) calls subroutines 
GRA and PPLOT for appropriate curves and drawings. 

Subroutine SUBR. - This subroutine is called by RUNGEK to calculate values of the 
derivatives for equations (2) to (5). 

Subroutine PPLOT . - This subroutine prepares data for shaft drawings. 

Subroutine GRA . - This subroutine prepares data for the following curves: 

(1) Deflection 

(2) Bearing against critical speed 

Subroutine RUNGEK . - This subroutine (1) performs the numerical integration and 
(2) saves numerical values of Y at every step for cases 1 and 2. 
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Subroutine CALPLT . - This subroutine plots the shaft and the superimposed deflec- 
tion curves and calls CALTIT. 

Subroutine CALTIT. - This subroutine writes Hollerith information and numbers on 
the plots. 


FORTRAN Symbols for NAMELIST Card 

The FORTRAN symbols for the NAMELIST cards are listed and defined in this sec- 
tion. The underlined symbols are those initialized by the program. 

A array for mass moment of inertia of disk about its axis of symmetry 

AI array for moments of inertia of sections 

AREA array for areas of sections 

BB array for mass moment of inertia of disk about axis through center of gravity 

and normal to axis of symmetry 

BRG array for bearing spring constants: 

BRG=0 - rigid bearing 
BRG=- 1 - coupling 

DEN array for mass density of disks 

DIA array for inside diameter of disks 

DRPM set equal to 200 

DX array for lengths of shaft sections 

E array for modulus of elasticity of shaft sections 

IC switch for boundary conditions at x = 0: 

IC = 1 - free overhang 
IC=2 - fixed end 

ICM switch used for plotting: 

ICM-Q - plot desired 
ICM=1 - no plot 

ID array for inside diameters of shaft sections 
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IH,KA 


IS 

ISK 


JXZ 


KA 

KB 

KL 

L 

LOC 

MOD 


ND 

NNN 

OD 

ODIA 

RHO 

RPM 


switches used for reading and writing data: 

IH=1 and KA=0 - data to be read in with FORMAT 300 
IH=1 and KA=1 - data to be read in with FORMAT 400 
IH=0 and KA=0 - data already read in with FORMAT 300 or to be entered 
through NAMELIST card 

IH=0 and KA=1 - data already read in with FORMAT 400 or to be entered 
through NAMELIST card 

used as a subscript for BRG array; tells which bearing will be used to form 
the x-axis of the bearing -stiffness -against -critical -speed plot; IS=1 

switch for the gyroscopic effect: 

ISK=0 - no gyroscopic data 

ISK=1 - gyroscopic effect to be considered 

ISK=2 - gyroscopic effect to be considered; data in computer 

switch used to write out data for the rpm -against -excess -moment curve: 
IXZ=0 - no writeout 
IXZ=1 - data desired 

see IH 

total number of bearings and couplings 

array for switches that tells program how gyroscopic data will enter: 

KL=4) - THIC, DEN, DIA, and ODIA 
KL-1 - A and BB 

number of shaft sections 

array for locations of disks for gyroscopic effect 

switch used to modify shaft data: 

MOD=0 - no modifications 

MOD=l - modifications to be made; program reads in shaft DATA cards 
and returns to beginning of program 

number of disks 

used to identify run number for critical -speed -against -bearing -stiffness plot; 
NNN=1 

array for outside diameters of shaft sections 
array for outside diameters of disks 
array for mass densities of shaft sections 
starting rpm (revolutions per minute); RPM= 1 
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upper limit of rpm interval; RPMF=50000 

percent used to determine maximum step size in integration; 
maximum step - length of shaft multiplied by RU ; RU-- Of 

variable used for scaling plots: 

SCALE=1 - full scale 
SCALE=2 - double size 
SCALE=. 5 - one-half scale, etc. 

(Because of the limitations of the Calcomp plotter, a maximum of a 5 -inch 
radius can be drawn. ) 

array for extra masses 

array for station numbers of bearings and couplings 
array for inside diameters of tapered sections 
array for outside diameters of tapered sections 
array for thicknesses of disks 

DATA INPUT 

All computer runs require a title card, a NAMELIST card, and a set of shaft data 
cards. The title card is used for problem identification on the printed and plotted output. 
The NAMELIST card contains all the input parameters and program switches. The shaft 
data cards contain the physical and geometric properties of the shaft. 

Two formats, 300 and 400, are used to read in the shaft data cards. FORMAT 300 
is used for the outside and inside diameters of each shaft section, along with the physical 
properties of the shaft. FORMAT 400 is used for the cross-sectional areas and the mo- 
ments of inertia of each shaft section, along with the physical properties of the shaft. 

See the section Preparation of Shaft Data Cards. 

The order in which the title card, NAMELIST card, and shaft data cards enter the 
computer is illustrated by figure 1. 


RPMF 

RU 

SCALE 


ss 

STA 

TAID 

TAOD 

THIC 
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Computer Output 


The program produces printed and plotted output. The printed output gives all the 
input data and the critical speeds. The program will also print the data necessary to plot 
the rpm -against -excess -moment curve. This option is used to confirm the converged 
critical speeds and is automatically printed whenever a discontinuity occurs in the excess 

moment. 

The plotted output is not only beneficial for analysis but also for simplicity in record- 
keeping. The plotting routine produces a drawing of the shaft with superimposed deflec- 
tion curves at the critical speed together with all pertinent information related to the 
shaft. Another form of plotted output is the critical -speed -against -bearing -stiff ness plot. 
For illustrations of the output, see the section Sample Problems with Output. 

In addition to the output already described, the program also prints the following 

error messages: ^ 

(1) Inside diameter is greater than outside diameter - shaft data card * . 

(2) Outside diameter of disk at station * is equal to or less than zero. 

(3) Location of disk is out of range. 

(4) Station location is out of range. 

(5) No convergence on critical speed occurred after 35 iterations. 


Preparation of Shaft Data Cards 

The shaft cards contain the physical and geometric properties of the shaft. The pro- 
cedure for the preparation of these cards for both tapered and untapered shafts follows: 
Shaft without taper. - Figure 2 is a scaled drawing of a multispan shaft system of a 
two -stage turbine. The system contains four bearings and two couplings. Also, there 
are eight extra masses^ (nonshaft mass) and one disk. To prepare this shaft system for 
the computer program, the following procedure is used: 

Step 1: Vertical lines are drawn perpendicular to the shaft where bearings, cou- 
plings, extra masses, and disks are located. 


Asterisk denotes station number or card number to be printed by program. 

^The program does not allow for an extra mass, bearing, or disk that is located at 
the initial end of the shaft x = 0. To achieve the equivalent, simply make section 
0-to-l small (0. 00001) and place the bearing, disk, or extra mass at station 1. 
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[Mil «■ Nonshaft mass 


Figure 2. - Multispan shaft system for a two-stage turbine (nontapered shaft). 


Step 2: Vertical lines are drawn perpendicular to the shaft where a change in every 
cross section of the shaft occurs, including the points x = 0 and x = L. 

Step 3: Starting at the left end of the shaft, number the vertical lines 0 to n. The 
numbered vertical lines represent stations, and the segments between stations are sec- 
tions. (Figure 3 illustrates this procedure for the first 25 sections of the shaft system 
of figure 2. ) 

Step 4: Fill out the shaft data sheet (fig. 4) with the appropriate values for each sec- 
tion and station. 

Shaft with taper. - Figure 5 is a scaled drawing of a three-span shaft system of a 
seven-stage compressor with a tapered portion. The system contains two bearings and 
15 extra masses. 

For the tapered portion the outside and inside diameters of the left end are placed in 
columns 6 to 10 and 11 to 15, respectively, of the shaft data sheet (fig. 4). The outside 
and inside diameters of the right end are placed in columns 26 to 30 and 31 to 35, respec- 
tively. Figure 6 illustrates the procedure for the first 10 sections of the shaft system of 
figure 5, and figure 7 is the data sheet. 

Blade preparation . - Because of the many variations in blade configurations, no set 
procedure is presented for the breakdown of a blade for proper input to the computer 
program. 

In general, a blade is broken down into different sections; and each section contains 
the cross-sectional area, moment of inertia, modulus of elasticity, and mass density. 
Figure 8 shows a breakdown of a blade, and figure 9 is the blade data sheet. 
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Bearing 




Shaft 

section 

Change in 
axial 
position, 
x, 
in. 

Outside 

diameter, 

in. 

Inside 

diameter, 

in. 

Modulus 

of 

elasticity, 
Ib/in. 2 

Mass 
density, 
Ib-sec'W 4 

Tapered 

outside 

diameter, 

in. 

Tapered 

inside 

diameter, 

in. 

Extra mass, 
lb-sec 2 /in. 2 
at station- 

Station 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

0-1 



1 




1 

2 

5 







2 

8 

0 

+ 

6 

7 

3 

3 

- 

5 
















1 

1-2 



7 

5 



1 

5 

2 

5 







1 


+ 






- 

















2 

2-3 


1 

5 




1 

5 












+ 






- 

















3 

3-4 



0 

5 



2 

1 












+ 






- 

















4 

4-5 


1 

3 




5 

1 












+ 























5 

5-6 


1 

2 




3 

3 












+ 






- 

















6 

6-7 



0 

8 



2 

9 












+ 






- 

















7 

7-8 



2 




1 

7 












+ 






- 

















8 

8-9 


1 

2 

5 



1 

5 












+ 






- 

















9 

9-10 


1 

2 

7 



1 

5 












+ 






- 

















10 

10-11 



1 




1 

5 





3 







+ 






- 

















11 

11-12 



6 

3 



1 

9 





4 

5 






+ 






- 












0 

0 

5 

0 

6 

12 

12-13 



3 

8 



1 

9 





4 

5 






+ 






- 

















13 

13-14 



0 

5 



1 

9 





6 

8 






+ 






- 

















14 

14-15 



2 




1 

9 




1 

3 

8 






+ 






- 

















15 

15-16 



1 

5 



1 

9 





4 

5 






+ 






- 

















16 

16-17 



1 

5 



1 

4 





4 

5 






+ 






- 

















17 

17-18 


2 

1 




3 

3 





4 

5 






+ 






- 

















18 

18-19 



5 




3 

3 





6 







+ 






- 

















19 

19-20 



1 




3 

3 




1 








+ 






- 

















20 

20-21 



3 

5 



3 

3 












+ 






- 

















21 

21-22 



9 




1 

5 












+ 






- 

















22 

22-23 



9 




1 

2 












+ 






- 

















23 

23-24 



3 




3 

1 

1 











+ 






- 

















24 

24-25 


1 





3 

1 

1 











+ 






- 

















25 

25-26 



















+ 





- 

















26 

26-27 



















+ 





- 

















27 

27-28 



















+ 





- 
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Figure 4. - Shaft data sheet for two-stage turbine (nontapered shaft). 






Extra mass 



Figure 6. - Detailed breakdown of first 10 sections of seven-stage compressor with tapered portion (fig. 5). 


Extra mass 




Shaft 

section 

Change in 
axial 
position, 
X, 
in. 

Outside 

diameter, 

in. 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

0-1 



7 

5 



2 

3 

8 



1-2 


3 

5 

3 



3 





2-3 


2 

5 

9 



3 





3-4 


1 

5 




4 

3 

6 



4-5 



4 




7 

2 

5 



5-6 



1 

5 



7 

2 

5 



6-7 



5 

3 



9 

4 




7-8 



2 

9 



9 

4 




8-9 



9 

4 


1 

0 





9-10 


5 

0 

3 


1 

0 

5 




10-11 












11-12 












12-13 












13-14 












14-15 












15-16 












16-17 












17-18 












18-19 












19-20 












20-21 












21-22 












22-23 












23-24 












24-25 












25-26 












26-27 












27-28 















Figure 7. - Shaft data sheet for seven-stage compressor with tapered portion. 





































































Sample Problems with Output 


The purpose of this section is twofold: (1) to show the flexibility of the program, and 

(2) to illustrate specifically the preparation of the NAMELIST card for particular prob- 
lems. The following six problems will be considered: 

(1) Multispan shaft with bearings and couplings, as well as a disk for the gyroscopic 

effect, and extra masses 

(2) Successive cases using same shaft data but different bearing constants 

(3) Problem 1 with modified shaft data 

(4) Free -free, single -span shaft 

(5) Set of cases to produce bearing-stiffness-against-critical-speed plot 

(6) Blade calculation with fixed-end boundary condition (cantilever beam and non- 
rotating uncoupled bending vibration) 

The program has built-in standard options; and consequently, they need not be 
placed on the NAMELIST card when a specific case is run. These standard options are 
indicated by an underscore in the section FORTRAN Symbols for NAMELIST Card. For 
example, ICM=0 is underscored and therefore indicates that a plot will always be part of 
the output. If, however, no plot is desired, ICM=1 must be put on the NAMELIST card 
for the specific problem being run. 

Problem 1 - multispan shaft. - The NAMELIST card described here contains the 
proper variables and constantsthat will enable the program to compute a critical speed 
for the shaft system in figure 2. The procedure for problem 1 is as follows: 

Title card TWO -STAGE TURBINE 

NAMELIST card $ SHAFT L=51, KB=6, SCALE=.25, STA=2, 9, 24, 34,41,45, 

BRG=2*50. E+03, 2* -1, 500. E+03, 100. E+03, ISK=1, ND=1, 
KL=1, A=. 80994, BB=. 79379, LOC=50$ 

Shaft data cards obtained from shaft system (fig. 2) 


where 

L=51 

KB=6 


SCALE=. 25 

ST A-2 , 9 , 24 , 34 , 4 1 , 4 5 


number of shaft sections 

total number of bearings and couplings 

number used for scaling plots 

station numbers of bearings and couplings 
STA(1)=2 
STA(2)=9 
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BRG=2* 50. E+03, 2* -1, array for bearing spring constants 

500. E+03, 100. E+03 BRG(1)=50. E+03 denotes flexible bearing at station 2 

BRG(2)=50. E+03 denotes flexible bearing at station 9 
BRG(3)=-1 denotes coupling at station 24 
BRG(4)=-1 denotes coupling at station 34 
BRG(5)=500. E+03 denotes flexible bearing at station 41 
BRG(6)=500. E+03 denotes flexible bearing at station 45 


ISK=1 


ND-1 


KL=1 


A=. 


BB= 


gyroscopic effect to be considered 
one disk for gyroscopic effect 

switch that tells that data for gyroscopic effect will enter with 
values of A and BB known 

array for mass moment of inertia of disk about its axis of 
symmetry 

array for mass moment of inertia of disk about axis through 
center of gravity and normal to axis of symmetry 


LOC=50 


station number where disk is located 


The following variables and constants need not be placed on the NAMELIST card because 
they are set by the computer program: DRPM, ICM, IC, RPM, RPMF, IH, and KA. If 
the program is to compute the values of A and BB, the NAMELIST card should be 

$ SHAFT L=51, KB =6, SCALE=.25, STA=2, 9, 24, 34,41,45, 

BRG=2* 5. E+04, 2* -1, 5. E+05, 1. E+05, ISK=1, ND=1, LOC=50, 

ODIA=5. 5, DEN=. 266E03, DIA=0, THIC = 1.5$ 

where ODIA, DIA, THIC, and DEN are the outside and inside diameters, the thickness, 
and the mass density of the disk. 

Figures 10 and 11 are the printed and computer -plotted output of problem 1. Fig- 
ure 12 is a hand-drawn plot of excess moment as a function of rpm, for problem 1. 
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TWO- stage turbine 


X 

OD 

ID 

EW 

e 

RHO 

TAOD 


C.1000 

1.2500 

-0. 

-0. 

0.28000E+08 

0. 73300E-03 

-0. 

-0. 

C.7500 

1.5250 

-0. 

-0. 

0.28000E+08 

0. 73300 E— 03 

-0. 

-0. 

1.5000 

1.5000 

-0. 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

C.0500 

2.1000 

-0. 

-0. 

0.28000E+08 

O.73300E-03 

-0. 

-0. 

1.3000 

5.1000 

-0. 

-0. 

0.28000E+08 

0. 73300E-0 3 

-0. 

-0. 

1.2000 

3.3000 

-0. 

-0. 

0.26000E+08 

0.73300E-03 

-0. 

-0. 

0.0800 

2.9000 

-0. 

-0. 

0.28000E+08 

0.73300E— 03 

-0. 

-0. 

C.2000 

1.7000 

-0. 

-0. 

0.28000E+08 

0.73300E— 03 

-0. 

-0. 

1.2500 

1.5000 

-0. 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

1.2700 

1.5000 

-0. 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

0.1000 

1.5000 

0.3000 

-0. 

0.28000E +08 

0.73300E— 03 

-0. 

-0. 

0.6300 

1.9000 

0.4500 

0.005060 

0.28000E+08 

O.73300E-03 

-0. 

-0. 

0.3800 

1.9000 

0.4500 

-0. 

0.28000E+08 

0.73300E— 03 

-0. 

-0. 

0.0500 

1.9000 

0.6800 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

0.2000 

1.9000 

1.3800 

-0. 

0.28000E+08 

0.733006—03 

-0. 

-0. 

0.1500 

1.9000 

0.4500 

-0. 

0.10000E+08 

0.25900E-03 

-0. 

-0. 

0.1500 

1.4000 

0.4500 

-0. 

0.10000E+08 

0.25900E-03 

-0. 

-0. 

2.1000 

3.3000 

0.4500 

-0. 

0. 10000E+08 

0.25900E-03 

-0. 

-0. 

0.5000 

3.3000 

0.6000 

-0. 

0.10000E+08 

0.25900E— 03 

-0. 

-0. 

0.1000 

3.3000 

1.0000 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

0.3500 

3.3000 

-0. 

-0. 

0.28000E+08 

0.733006-03 

-0. 

-0. 

0.9000 

1.5000 

-0. 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

0.9000 

1.2000 

-0. 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

0.3000 

3.1100 

-0. 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

1.0000 

3.1100 

-0. 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

0.6500 

4.7500 

2.2200 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

0.5100 

3.1800 

1.3600 

-0. 

0. 280006+08 

0.733006-03 

-0. 

-0. 

0.7900 

2.4000 

1.3600 

-0. 

0.28000E+08 

0.73300E— 03 

-0. 

-0. 

14.5000 

2.0000 

1.3600 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

0.1700 

2.0000 

1.3600 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

0.5200 

3.0000 

1.3600 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

0.7400 

3.1900 

1.3600 

-0. 

0.28000E+08 

0.73300E— 03 

-0. 

-0. 

0.6000 

4.7500 

2.2500 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

0.4900 

3.1000 

0.8800 

-0. 

0.28000E+08 

0.733006-03 

-0. 

-0. 

0.7500 

3.1000 

0.8800 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

1.0200 

1.5600 

0.8800 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

0.6000 

2.0000 

0.8800 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

0.2800 

2.1500 

0.8800 

0.001400 

0.Z8000E+08 

0.733006-03 

-0. 

-0. 

C.2800 

2.1500 

0.8800 

-0. 

0.28000E+08 

0. 733006—03 

-0. 

-0. 

1.1800 

2.5000 

0.8800 

0.000643 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

C.6600 

2.5000 

0.8800 

-0. 

0. 28000E+08 

0.733006-03 

-0. 

-0. 

4.1000 

2.5000 

0.8800 

0.000908 

0.28000E+08 

0.73300E— 03 

-0. 

-0. 

2.6800 

2.5000 

0.8800 

-0. 

0.28000E+08 

0. 733006-03 

-0. 

-0. 

C.860G 

2.5000 

0.5000 

-0. 

0.28000E+08 

0.73300E— 03 

-0. 

-0. 

0.6300 

2.5000 

1.3100 

-0. 

0.28000E+08 

O.73300E— 03 

-0. 

-0. 

C.2000 

2.5000 

1.3100 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

0.6600 

2.5000 

1.6600 

0.000643 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

1.2500 

2.5000 

1.6600 

0.004300 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

0.4000 

2.5000 

-0. 

0.015000 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

1.1400 

2.5000 

-0. 

0.013700 

0.28000E+08 

0.73300E-03 

-0. 

-0. 

2.2200 

2.5000 

-0. 

-0. 

0.28000E+08 

0.73300E-03 

-0. 

-0. 


BRG AT STATION 2K- 
BRG AT STATION 9K* 
COUPLING AT STATION 
COUPLING AT STATION 
BRG AT STATION 41K- 
BRG AT STATION 45K- 


BOUNDARV CONDITIONS- 
0.500000E+05 
0.500000E+05 
24 
34 

O* 500000E+06 
0. 100000E+06 


•SHEAR 


ANO MOMENT EQUAL 


ZERO AT X-0 AND 


TOTAL MASS • 0. 19590678E+00 

THIS RUN CONSIDERS GYROSCOPIC EFFECT 


DISK AT STATION 50 A-B- 0.16I500E-01 


A DISCONTINUITY OCCURS AT 15000. IN THE EXCESS MOMENT VS RPM CURVE 

DATA FOR THIS CURVE FOLLOWS 


1ST CRITICAL SPEED - 
2ND CRITICAL SPEED - 


0. 0.42578142E+16 
SCO. C.63256416E+08 
1600. 0.13462204E+08 
24C0. 0.43670 5 TOE +07 
32C0. 0. 13291 924E *07 
40C*. 0.75740273E+05 
48C0.-0. 45248 193E+06 
56C0 .*0.62327 2 3 IE +06 
64CC.-0.5950C831E+06 
72C O.-O. 44778 150E +06 
8CCO.-0.22T71575E+06 
8800. 0.35T67375E+05 
96C0. 0. 3 2 26 7 92 5E +06 
104C0. 0.61957 6 62E +06 
112CC. 0.91 9364 87E +06 
120C0. 0. 12245032E+0T 
128C0. 0.1559542 3E +07 
136C0. 0.20352 74 7E +07 
144C0. 0.34942161E+07 
152C0.-0.257937T6E+07 


200. 0.1060107BE+10 
1000. 0. 39330808E+08 
1800. 0.100D2868E+08 
2600. 0.33179008E+07 
3400. 0.91103281E+06 
4200.— 0. 1040759 5E+06 
5000.-0. 52044806E+06 
5800.-0. 6 30 8936 2E +06 
6600.— 0.5 671 5800 E+06 
7400.-0. 39823638E+06 
8200.-0. 16500613E+06 
9000. 0. 10592550E+06 
9800. 0.39636200E+06 
10600. 0.69436575E+06 
11400. 0.99473450E+06 
12200. 0. 13036270E+07 
13000. 0.16570448E+07 
13800. 0.22210095E+07 
14600. 0.49467520E+07 
15400. -0.98175519E+06 


400. 0.26270325E+09 
1200. 0. 26344510E+08 
2000. 0. 754058 60E+07 
2800. 0.24977620E+07 
3600. 0.57266344E+06 
4400.— 0.24 8433 30E +06 
5200.-0. 56997019E+06 
6000.-0. 62790281E+06 
6800.-0. 53286000E+06 
7600. -0.344737386+06 
0400.— 0. 99979 125E +05 
9200. 0. 17727375E+06 
10000. 0. 470487 12E+O0 
10800. 0.76925337E+06 
11600. 0. 107051 86E + 07 
12400. 0.13851142E+07 
13200. 0.17646734E+07 
14000. 0.247234126+07 
14800. 0. 11623892E+08 
15600. -0.31028778E+06 


3RD CRITICAL SPEED - 
600. 0.114965116+09 
1400. 0.18525405E+08 
2200. 0.57310430E+07 
3000. 0. 1 8484066E+07 
3800. 0.29823668E+06 
4600.-0. 36304434E+06 
5400.-0. 603549066 +06 
6200.-0. 615575006+06 
7000.-0. 492850 13E+06 
7800.-0. 28775487E+06 
8600.-0. 3Z956375E +05 
9400. 0. 249591 13E +06 
10200. 0. 54492625E+06 
11000. 0. 84423800E +06 
11800. 0.1 1469808E +07 
12600. 0. 147004176+07 
13400. 0.18879704E+07 
14200. 0. 28471443E+07 
15000.-0. 11240952E+08 
15800. 0. 57058683E+05 


4084. 

8696. 

15769. 


Figure 10. - Printed output for problem 1 (multispan shaft with bearings and couplings, as wen as a disk for the gyroscopic 
effect, and extra masses). 



•01 ajnBij eas ‘pajja oidoosojAB ;o s|ie}9p joj -906561 '0 ‘sseui |eioi ’i Luajqojd joj indjno papoid - ‘II 9jn6y 



DEFLECTION. 

-2D. 00 -16.00 -12.00 -8. DO ->i.00 0.00 H.00 

sK — * 1 ■ x — m RING- 50000. 


.6 


5x: 



Speed, rpm 

Figure 12. - Hand-drawn plot of excess moment as function of rpm, for problem 1. 

Problem 2 - successive cases using same shaft data . - The procedure for problem 2 
is as follows: 

For case 1: 

Title card same as for problem 1 

NAMELIST card same as for problem 1 

Shaft data cards same as for problem 1 

For case 2: 

Title card TWO-STAGE TURBINE BRG AT STA2=6. E+06 

NAMELIST card $ SHAFT IH=0, SCALE=.25, ISK=2, BRG(1)=6. E+06$ 

where 

IH=0 switch that indicates shaft data are in computer 

ISK=2 switch that indicates gyroscopic data are in computer 
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BRG(1)=6. E+06 


new spring constant at station 2 (all other BRG values remain 
unchanged) 

Problem 3 - problem 1 with modified shaft data . - A change from 2. 5 inches to 
3. 0 inches will be made on the outside diameter for sections 39-40 and 50-51 inclusive. 
The procedure for problem 3 is as follows: 

Title card TWO-STAGE TURBINE - ORIGINAL DATA 

NAMELIST card $ SHAFT L=51, MOD=l$ 

Shaft data cards Same as for problem 1, where MOD=l is the switch that re- 

turns flow of program to beginning, after shaft data cards 
are read 

TWO -STAGE TURBINE - MODIFIED SHAFT DATA 

Same as for problem 1, plus IH=0, OD(40)=12* 3. $ 

outside diameter of 3 inches for 12 consecutive sections start- 
ing at section 39-40 

Problem 4 - free -free, single -span shaft . - The procedure for problem 4 is as 
follows: 

Title card TWO STAGE TURBINE - NO BEARINGS OR COUPLINGS 

NAMELIST card $ SHAFT L=51, KB=0, SCALE=.25$ 

Shaft data cards Same as for problem 1 

Problem 5 - bearing -stiffness -against -critical -speed plot . - The NAMELIST cards 
for this problem illustrate the ability of the program to produce a bearing-stiffness - 
against -critical -speed plot. This plot requires six computer runs that use the same 
shaft data cards. The scale on the x-axis of this plot is determined by changing the 
spring constant of the same bearing for each run. The first computer run has the 
smallest bearing constant, the second has the next larger, and so forth. The data used 
for these runs represent the shaft system of figure 5. 

Case 1: 

Title card SEVEN-STAGE COMPRESSOR - CASE 1 

NAMELIST card $ SHAFT L=39, STA=2,31, BRG=2* 300. E+03, SCALE=.25, 

KB=2$ 

Shaft data cards obtained from shaft system of figure 2 


Title card 
NAMELIST card 
OD(40)=12* 3. 
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CASE 2 

$ SHAFT IH=0, NNN=2, BRG=2* 500. E+03 , SCALE-. 25$ 


Case 3: 

Title card 
NAMELIST card 

Case 4: 

Title card 
NAMELIST card 

Case 5: 

Title card 
NAMELIST card 

Case 6: 

Title card CASE 6 

NAMELIST card $SHAFT IH=0, NNN=6, BRG=2* 200. E+04 , SCALE=.25$ 

where NNN is the number used to identify run number and is used only when a bearing- 
stiff ness -against -critical -speed plot is desired. Figures 13 and 14 are printed and 
computer -plotted output of case 1. Figure 15 is the bearing-stiffness -against -critical - 
speed plot for cases 1 to 6. 

Problem 6 - blade calculation with fixed end boundary condition . - This problem in- 
volves a cantilever beam having nonrotating, uncoupled bending vibration. The procedure 
for problem 6 is as follows: 

Title card BLADE CALCULATION - FIXED END 

NAMELIST card $SHAFT L=20, KB=0, IC=2, KA=1, DRPM=1000. , RPMF=75000. $ 

where 

IC=2 switch for fixed end boundary condition at x = 0 

KA=1 data read in by FORMAT 400 

Blade data cards obtained from blade illustrated in figure 8 

Figure 16 is the printed output of problem 6. 


CASE 5 

$ SHAFT IH=0, NNN=5, BRG=2* 150. E+04 , SCALE=.25$ 


CASE 4 

$ SHAFT IH=0, NNN =4, BRG=2* 100. E+04 , SCALE= 25$ 


CASE 3 

$SHAFT IH=0, NNN=3, BRG=2*700. E+03, SCALE=. 25$ 


Case 2: 

Title card 
NAMELIST card 
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TEST 

PRGELEN 

TAPERED SHAFT 

X 

DD 

ID 

EW 

C.75CC 

2.3800 

1.3800 

0.014800 

3. 53CC 

3.0000 

0.7500 

-0. 

2.59CC 

3 .0000 

0.7500 

-0. 

1.5CCC 

A. 3600 

0.7500 

-0. 

C.4CCC 

7.2500 

0. 7500 

-0. 

C. 15C0 

7.2500 

0. 

-0. 

C .5300 

5.4000 

0. 

-0. 

C.29CC 

9.4000 

1.8600 

0.015500 

C.94CC 

1C. 0000 

8.6200 

0.076800 

5.0300 

1C. 5000 

9.8800 

0.113900 

4.65CC 

12.3800 

11.8000 

0.116800 

3.8600 

14.0000 

13.5000 

0.103400 

3.0CCC 

15.2000 

14.7000 

0.086400 

2.50CC 

16.2000 

15.7000 

0.086600 

2.17CC 

16 . 8400 

16.2400 

0.002500 

C.25CC 

16.8400 

9.2400 

-0. 

C.410C 

15.9600 

9.2400 

0.015200 

C.30CC 

12.5000 

9.2400 

0.105200 

C . 7800 

11.1600 

9.2400 

-0. 

C.27C0 

11.1600 

8.1600 

-0. 

C • 2400 

11.1600 

7.2200 

-0. 

C.07C0 

11.1600 

6.1000 

-0. 

C . 1500 

5.1200 

6.1000 

-0. 

C.06CC 

6. 4400 

6.1000 

-0. 

C . 2800 

8.4400 

3.6200 

-0. 

C • 2400 

7.9000 

3.6000 

-0. 

1 • 2700 

7.5000 

3.6200 

-0. 

C . 5800 

5.5600 

3.6200 

-0. 

l.OOCC 

5.5600 

2.2500 

0.014400 

1.1C00 

4.0000 

2.2500 

-0. 

1.84CC 

4.0000 

-0. 

-0. 

1.51CC 

4.0000 

-0. 

-0. 

1.9900 

3.9200 

-0. 

-0. 

1.1300 

3.2500 

-0. 

0.055400 

C • 8400 

3.2500 

-0. 

-0. 

C. 9600 

3.0000 

-0. 

0.010200 

1.28CC 

3.0000 

-o. 

0.006700 

C.81C0 

3.0000 

-0. 

0.079200 

C.76CC 

3.0000 

-0. 

-0. 


E 

0.29000E+08 
0.29000E+08 
0.29000E+08 
0.29000E+08 
0.29000E+08 
0 .290006 + 08 
0.29000E+08 
0.2 9000E+08 
0. 29000E+08 
0.29000E+08 
0.29000E+08 
0.29000E+08 
0.29000E+08 
0.29000E+08 
0.29000E+08 
0.29000E+08 
0.29000E+Q8 
0.29000E+08 
0 • 29000E+08 
0.29000E+08 
0.29000E+08 
0.29000E+08 
0.29000E+08 
0.29000E+08 
0.29000E+08 
0. 29000E -KJ8 
0.29000E+08 
0.29000E+08 
0.29000E+08 
0.29000E+08 
0.29000E+08 
0.29000E+08 
0.29000E+08 
0.29000E+08 
0.2 9000E+08 
0 . 29O0OE + 08 
0 • 29000E + 08 
0.29000E+08 
0.29000E+08 


RHO 

0.74000E-03 -0. 

0.74000E-03 -0. 

0.74000E-U3 -0. 

0.7 4000E-03 -0. 

0.74000E-03 -0. 

U.74000E-03 -0. 

0 • 74000E-0 3 -0. 

0.74000E-03 
G.74000E-03 
0.74000E-03 
0.74000E-03 
0. 74000E-0 3 
0.74000E-03 
O.74000E-03 
0.7 4000E— 0 3 
0.74000E-03 -0. 

0 • 74000E-G3 -0. 

0.74000E-03 -0. 

0. 74000E-03 -0. 

0. 74000E-0 3 -0. 

G.74000E-03 -0. 

0.74000E-03 -0. 

0.74000E-03 -0. 

0.74000E-03 -0. 

O.74000E— 03 -0. 

0. 74000E-03 -0. 

0. 74000E— 03 -0. 

0. 74000E-03 -0. 

0.74000E-03 -0. 

0.74000E-03 -0. 

0. 74000E— 03 -0. 

0. 74000E-03 -0. 

J.74000E-03 -0. 

0.74000E-03 -0. 

0.7 4000E-03 -0. 

0.74000E-03 -0. 

0.74000E-03 -0. 

0.74000E-03 -0. 

0. 74000E-03 -0. 


TAI7D TAID 

- 0 . 

- 0 . 

- 0 . 

- 0 . 

- 0 . 

- 0 . 

- 0 . 

- 0 . 

9.8800 
11.8000 
13.5000 

14.7000 

15.7000 
16.2400 
- 0 . 

- 0 . 

- 0 . 

- 0 . 

- 0 . 

- 0 . 

- 0 . 

- 0 . 

-0. 

- 0 . 

-0. 

-0. 

-0. 

-0. 

-0. 

-0. 

-0. 

-0. 

- 0 . 

-0. 

-0. 

-0. 

- 0 . 

- 0 . 

- 0 . 


- 0 . 

10.5000 

12.3800 

14.0000 

15.2000 

16.2000 
16.8400 
- 0 . 


BOUNDARY CONDITIONS SHEAR AND MOMENT EQUAL ZERO AT X=0 AND 

BRG AT STATION 2K= 0.300000E+06 

BRG AT STATION 31K= 0.300000E+06 

TOTAL MASS = 0. 14907389E+01 


1ST CRITICAL SPEED = 
2ND CRITICAL SPEED = 
3RD CRITICAL SPEED = 


5694. 

9031. 

25432. 


Figure 13. - Printed output for problem 5 (bearing-stiffness-against-critical-speed plot), case 1. 




Figure 14 - Plotted output for problem 5, case L Total mass, 1.490739. 


TfiH MfiSSn .055400 




TEST PROBLEM TAPERED SHAFT 



Q FIRST CRITICAL SPEEO 
A SECOND CRITICAL SPEED 
+ THIRD CRITICAL 5PEED 
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8D 120 

BEARING STIFFNESS 


Figure 15. - Bearing stiffness as function of critical speed for all six cases. 


BLACE CALCULATION CANTILEVER BEAM — FIXED END 


X 

AREA 

I 

EW 

E 

RHO 

0.06900 

0.490E-01 

0.178E-04 

-0. 

0.275006*08 

y. 750006-03 

C .06900 

0.477E-01 

0.161E-04 

-0. 

0.27500E*08 

0.75000E— 03 

C .06900 

0.465E-01 

0.147E-04 

-0. 

0.27500E+08 

0.75000E-03 

C .06900 

0.452E-01 

0.134E-04 

-0. 

0.27500E+08 

0.75000E - 03 

C. 06900 

0.440E-01 

0. 123E-04 

-0. 

0 .275006*08 

0. 75000E— 03 

C. 06900 

0.4006-01 

0. 113E-04 

-0. 

0.275006*08 

0.75000E-03 

C. 06900 

0.417E-01 

0.104E-04 

-0. 

0.2750DE+08 

0.75000E-03 

C .06900 

0.404E-01 

0.951E-05 

-0. 

0.275006*08 

0.75000E-03 

C. 06900 

0.391E-01 

0. 864E-05 

-0. 

0.27500E+08 

0.75000E-03 

C .06900 

0.378E-01 

0. 780E-05 

-0. 

0.275006*08 

0. 750006-03 

0.06900 

0.365E-01 

0.701E-05 

-0. 

0.275006*08 

0.75000E— 03 

C. 06900 

0.353E-01 

0. 630E-05 

-0. 

0.275006*08 

0.75000E— 03 

0.06900 

0.340E-01 

0. 564E-05 

-0. 

0. 27500E+0 8 

0.75000E— 03 

C. 06900 

0.327E-01 

0.503E-05 

-0. 

0.27500E+08 

0.75000E-03 

C. 06900 

0.3136-01 

0.446E-05 

-0. 

0.275006*08 

0.750006-03 

0.06900 

0.300E-01 

0.392E-05 

-0. 

0.27500E*08 

0. 75000E— 03 

0.06900 

0.286E-01 

0.342E-05 

-0. 

0.27500E+08 

0.750006-03 

0.06900 

0.272E-01 

0.296E-05 

-0. 

0. 275006*08 

0.75000E-03 

0.06900 

0.259E-01 

0.259E-05 

-0. 

0.275006*08 

0 • 75000E-03 

0.06900 

0.248E-01 

0.232E-05 -0. 0.27500E+08 

BOUNDARY CONDITIONS DEFLECTION 

0.750006-03 

AND SLOPE EQUAL ZERO AT X»0 


AND SHEAR AND FOMENT EQUAL ZERO AT X-L 

THIS RUN HAS NO BEARINGS OR COUPLINGS 


1ST CRITICAL SPEED = 


Figure 16. - Printed output for problem 6 (blade calculation with fixed end boundary condition). 


70135. 
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CONCLUDING REMARKS 


A FORTRAN IV computer program using a modified method for calculating critical 
speeds of rotating shafts has been described. Its main features are flexibility, plotted 
computer output, and a minimum amount of data preparation. The program is easy to 
use and is a useful tool for design engineers concerned with rotating machinery. 

The program can be adapted to any computer system that uses FORTRAN IV, but a 
plotting routine must be supplied by the user . 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, May 18, 1973, 

501-24. 
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APPENDIX A 


SYMBOLS 


2 

A mass moment of inertia of disk about its axis of symmetry, lb -in. -sec 

B mass moment of inertia of disk about axis through center of gravity and 

2 

normal to axis of symmetry, lb -in. -sec 


D. 


D. 


D o 

E(x) 

F(x) 

G(x) 

h 

I(x) 

K 

L 

M(x) 

M(x + ) 

M(x") 

M nl’ M n2’ M n3 

S n 

V(x + ) 

V(x“) 

V nl’ V n2’ V n3 


Y(x) 

Y g (x) 

Y „l’ Y n2’ Y n3 

7)(x) 

o(x) 


diameter of disk, in. 

inside diameter of shaft section, in. 

outside diameter of shaft section, in. 

2 

modulus of elasticity at x, lb/in. 

shear coefficient for shear deflection equation 

2 

modulus of rigidity at x, lb/in. 
thickness of disk, in. 

. 4 

diametral moment of inertia at x, in. 
bearing stiffness factor (spring rate), lb/in. 
length of shaft, in. 

moment for a general cross section, lb -in. 

moment to right of a station, lb -in. 

moment to left of a station, lb -in. 

moment for solutions 1, 2, and 3 at station n 

nonshaft mass at station n, lb -sec /in. 

shear to right of a station, lb 

shear to left of a station, lb 

shear for solutions 1, 2, and 3 at station n, lb 

station or axial position on shaft 

deflection for a general cross section, in. 

shear deflection for a general cross section, in. 

deflection for solutions 1, 2, and 3 at station n, in. 

2 2 

mass per unit length of a general cross section, lb -sec /in. 
slope at a general cross section (nondimensional) 



slope for solutions 1, 2, and 3 at station n (nondimensional) 
Poisson’s ratio 

2 4 

mass density, lb-sec /in. 

speed of rotation or frequency of vibration, rad/sec 
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APPENDIX B 


LISTING OF FORTRAN PROGRAM 


4**4********** I******************************************************** 

CRITICAL SPEED PROGRAM — CSP. 

READS ANC WRITES ALL INPUT DAT A* SETS 
SWITCHES TO REGULATE THE FLOW OF THE 
PROGRAM. CALCULATES GYROSCOPIC DATA (IF ANY) 

4*4 4 4 *******4 4 ******************* 4*4* ********************************** 

REAL ID(2CO), MAS$(200) 

INTEGER STA(IO) 

COMMON IC,MASS,STA 

CQMMCN/CLPLOT/XPEN»YPEN,NX,NY,I PEN , XL ABEL ( 10 ) , YL ABEL ( 10 ) 

COMMON RAT(10),CS(3>, DX ( 200 ) , 00 (200 ) , SS t 200 ) , E ( 200 ) , RH0( 200 ) , 
1TA0C(200) ,TAID(200) , AREA ( 200 ) , A I (200) , DEL (200) ,BRG( 10),BC( 10), 
2ZZ(2CO) ,EI (200) ,GAM(200> *LOC (200 ) , THIC ( 200 ) t OEN ( 200 ) ,ODI A( 200) , 

3 0 I A ( 200 ) *CRIT (6,3) , K0{ 10 ) ,DY(58),Xt 500 ) , Y1 ( 500 ) , Y2 ( 500 ) *KL ( 200 ) , 

4 EM (400 ) * VP ( 10 ) , SL ( 200 ) , SLL ( 200 ) * ON (10),TITLE(12) 

CCMMCN PI, KEY, NNN,IC, SCALE, ICM, I SK , I JJ,KA, INN, IXY # RPM , NO,NUM, 

1TP ,RP , KB , WW,IR,J,KK,DRPM,FF,L , TOTAL , IX Z , NOO , RPMF , MOD 
DIMENSICM A(200) ,88(200) 

NAME L I ST/SHAFT/DX,OD ,ID,S$,E ,RHO,T AOD,TA ID, AREA ,AI,BRG,STA,KB, 

IRPM , CRPM , IH, IC , SC ALE ,ICM, ISK ,KA, KL ,L *ND,LOC, THIC, DEN, ODI A,D1A, A, 
2BB,NC0,0EL,IS ,NNN, I XZ ,RPMF , MCD,RU 
30C FCRMAT(3F5.3,2E5.1,2F5.3,F5.5) 

AOC F0RMAT(FE.5,2E9.5,2E5.1,F6.5) 

7 FORMAT(1FO, 30X ,60HBOUNDARY CONDITIONS DEFLECTION AND SLOPE EQUAL 

1 ZERO AT X*0/39H AND SHEAR AND MOMENT EQUAL ZERO AT X»L I 
IOC F0RMAT(1FG,5X,1HX,9X, 2H0D » 8X ,2H ID,9X,2HEW,11X, 1HE, 12X , 3HRH0, I0X, 
14HTACD,1CX,4HTAID) 

11 FORMAT (1> ,3F10.4,F12.6,2E14.5,2F10.4) 

401 FORM A T ( 1F0 ,5 X , 1HX ,9X ,4HAREA , 8X , 1HI ,9X , 2HEW , 1 IX , 1HE, 1 2X, 3HRH0 ) 

402 FORMAT ( 1 > , F 10. 5, 2E12. 3, F 12. 6 ,2E14. 5 ) 

40CC FORMAT { 1F0 ♦ 5X , 12HT0TAL MASS -E15.8) 

6 F0RMAT(1F1, 30X ,66H BOUND ARY CONDITIONS SHEAR AND MOMENT EQUAL 

1 ZERC AT X*0 AND L) 

202 FORMAT ( IMG , 10X , 14HBRG AT ST ATI ON 13 , 2HK*E15 .6 ) 

20C FORMAT(1FC,10X,20HRIGID BRG AT STATI0NI4 ) 

203 FGRMAT(1F0,10X,20HC0UPLING AT STATI0NI4) 

25CC FORMAT ( 1F0 » IOX »62H INSIDE DIAMETER GREATER THAN OUTSIDE 0 I AMETER-SH 
1 AFT DATA CARD 15) 

25C1 F0RMAT(1F0,20X,37HTHIS RUN HAS NO BEARINGS OR COUPLINGS) 

2502 FORMAT ( IF 0 , IOX ,35H0UT$I DE DIAMETER OF DISK AT STATION, 14, 1H*,F9. 5) 
25C5 FORMAT ( 1F0 , IOX, 29HL0CATI0N OF OISK OUT OF RANGE , 2X , 8HSTAT ION* 14 ) 
25C2 FORMAT ( 1F0 , IOX ,29HSTATI ON LOCATION OUT OF RANGE , 2X, 8HSTAT ION-I 4 ) 
3CCC FORMAT (1F0,15X,36HTHIS RUN CONSIDERS GYROSCOPIC EFFECT) 

30C2 FORMAT ( IF 0 ,15X,15HDI SK AT ST ATI 0NI4,5X ,4HA-B*E15.6 ) 

1CCC FGRM AT ( 12 A6 ) 

20CC FORM A T ( IF 1 , 12 A6 ) 

PI*3. 1415927 
CC*P I /32. 

5 BG 590 IV=1,10 
K0< I V ) =0. 

59C RAT ( I V ) »C « 
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DC 868 1*1 , 3 
686 CS(I)=G. 

DC 332 1 = 1 #200 
392 KL(U=0 
13*1 
RL*.C1 
1 1 Z = C 

RPRF “5Q0C0 • 

MOD = C 
RPP=.l 
DRPP = 200 • 

KE Y=C 
NNN = 1 
IH*1 
IC*1 

5CALE*1. 

ICR*C 

isk*c 

ijj=i 

KA *0 
INN* 1 

rxY*c 

I X z *o 

RE AD ( 5 » 1CCC ) TITLE 
WRITE 16, 2CC0) TITLE 
READ 15, SFAFT) 

ON (NNN ) *ERG < I S ) 

rPtlH.EQ. l.AND.KA.EQ.O) GO TO 67 
IF ( IH.EQ.O* AND.KA.EQ.O) GO TO 78 
IF< IH.EQ.l.AND.KA.EQ.l) GO TO 79 
GO 10 87 

67 REACt 5 t 3 CCMDX(l)fOD<n*IDtl>»E<I) *RHO ( I ) , TAOD( II , TA I D< I) t SS < I ) , 
11*1, L) 

WRITE(6*11) <DX (I ) ,CD(I ) * I D (I I » SS <1 1 1 E< 1 ) ,RHOt I ) * T AODC I ) » 

1 T A I C ( I ) ,1*1, LI 
DO 919 1*1, L 

IF ( I D ( I )«GT»OD<I )*OR*TAIDU I .GT*TAOD( I I I GO TO 991 
©0 TO 913 

391 WRtTE(6,25C0) I 

IIZ*1 

313 CONTINUE 

IE ( 1 1 Z *EC * 1) GO TO 5 
GO TC 90 

73 READ <5,4C0 ) (DXCI ) , AREA ( I ) , A I ( I ) , E ( I ) »RHC ( I ),S$<I>tl*ltL) 

WRITE(6,A02) (DX(I),AREA(I ) »AH I ) ,SS( l ) y E ( I ) ,RHO( I ) , 1*1, L) 

ICP*1 

30 IE < HOD* EC • 1IG0 TO 5 
nccmrpmf-rpmi/drp* 

NO = NOO 
TP*C • 

DO 55 1*1, L 
55 TP*CX(IM TP 
STEP “TP *RL 
DO 351 1*1. L 

TFISTEP •LE.DX(I)) GO 10 «06 
DEL(I)=1. 

GO TC 951 
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806 I5*CX< I ) /STEP +.999 
DEL ( I ) -1 5 
951 CONTINUE 
SUP*C. 

DO 950 l*l*L 
95C SUP*CEL(I)+SUM 
NUP*SU*+1* 

IF<NNN.GT*1> GO TO 7C0 
00 44 3 1*1 » 6 
DO 443 I\*l,3 
442 CRITU,IV)*0. 

70C RP=RPM 
LP*KB+1 

DO 8C5 J2*LP » 10 
8RG( JZ)*C. 

805 STA ( JZ )=C 

GO TC( 1C » 2C ) « IC 
10 BC ( 3 ) =0. 

BC 14 ) -0. 

BC 1 5 I *0. 

BC<6 )*.0CQ1 
BC(7)*0. 

BC(fi)*0. 

8C<9)**CC01 
BC UC ) =0 • 

WRITE<6»6) 

GO TC 132 
20 BC<3)*.CCG1 

BC ( A )*0* 

BC ( 5 ) *0* 

BCt 6 1*0* 

BC 1 7 ) *0* 

Bcm*.ocoi 

eci9)*o* 

8C1 10*0* 

WRITE<6,7 ) 

132 IFIKB.EQ.C) GO TO 41 
DO 66 1=1 , KB 

IP ( BRG (I ) • EQ* 0*1 GO TO 43 
IF ( 8RG ( l ) • EC.-l* ) GO TO 44 
WRITE(6,202) STAU),BRG(I) 

IF t STA ( I ) *GT* L .CR.STAU ).LE.O) GO TC 666 
GO TO 66 

43 RRITE16.20C) STA C I 1 

IF<$TAtI).GT*L *0R* STA ( I ) • LE .0 ) GO TO 666 
GO TC 66 

44 WR ITE ( 6 » 203 ) STAU) 

IF { STA ( I ) *GT* L *OR.STA( D.LE.O) GO TO 666 
GO TC 66 

666 WRITE(6»25C2) STAU) 

IIZ*1 

66 CONTINUE 

IFUIZ.EC* 1) GO TO 5 
GO TC 77 

41 WRITE<6,2501) 

77 IFUSK-2) 8 2 * 38 1 82 
82 DO 8C KZ* 1 1 200 
80 ZZ ( K Z ) =0 • 

38 IF (KA* EQ. C) GO TO 33 
DO 62 1*1 *L 



t* ASS< I )=AREA ( I ) *RHO ( I ) 
ei ( n*E( n*Ai ( i) 

62 <5AMI>=0. 

GO TG 61 

33 CO 25 1=1,1 

IFUAOOm.EQ. 0.) GC TO 24 
SL<I)*<TAOD(I l-OCini/OXUl 
GO TC 25 

24 SL ( l ) =0* 

25 CONTINUE 

DC 26 1=1, L 

IF ( TAOD ( 1 ) .EQ. 0.) GO TO 23 
SLL m = CTAIO(I)-ID( I H/DXU ) 

GO TC 26 
23 SLL<I)=0. 

26 CONTINUE 
TGTAL*0. 

DO 556 1*1, L 
TIL=PI*REC I I )*DXU)/3. 

IF ( T AOD ( 1 ) .EQ.O. ) GO TO 555 

VX*T IL* (CD i I ) **2/4. +00 < I )*TA0D< I 1/4.+TA0CI I )**2/4. ) +SS( I ) 
TGTAL=T01AL+VX-TIL*< 10(1 )**2/4.+ID( I )*TAID( I )/4.+TAID( 1 1**2/ 4. ) 
GO TG 556 

555 TOTAL* (OC ( I )**2-ID ( I )**2)*PI *DX< I ) *RHO( I )/4.+T0TAL +SS(I) 

556 CONTINUE 

HR I TE < 6,4CC0 I TOTAL 
61 IF< ISK.EC.C ) GO TO 88 
IFUSK.EC. 2) GO TO 72 
DO 71 I =1 ,ND 
NA*LCC ( I ) 

IFILCC m.GT.L.OR.LOCII).LE.O> GO TO 748 
IEIKLt I ) .EC. 1) GO TO 48 
IF ( OD I A I Il.GT.O. ) GO TO 894 
WR I T£ ( 6 , 2 503 ) LOC ( I) , OD I A ( I ) 

GO TC 5 

894 D=CC*THIC ( I )*DEN (I) *(001 Ad) **4- DIA(I)**4> 

CD=.5*D*(1.+4.*THIC(I ) **2/ ( 3.* ( ODI A ( I ) **2- 01 A (I )**2) ) ) 

ZZ ( N A 1 *D-C 0 
GO TC 71 

48 ZZ(NA)=A(I )-BB ( I ) 

GO TC 71 

748 WRITE<6,25C5) LOC (I) 

IIZ*1 

71 CONTINUE 

IF { I I Z.EC . 1 ) GO TO 5 

72 fcR I TE ( 6 , 3CC0 ) 

00 73 1 = 1, NC 
NA*LCC ( I ) 

ViR ITE ( 6,3002 ) LOCU),ZZ(NA) 

73 CONTINUE 
88 CALL MAIN 
GO TC 5 
ENC 
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*********************************************************************** 

MAIN SUBROUTINE 

SETS UP INITIAL CONDITIONS FOR 

NUMERICAL INTEGRATION ANO THEN CALLS SUBROUTINE 
RUNGEK. PERFORMS CALCULATIONS FOR CHANGES 
BROUGHT ABOUT BY EXTRA MASSES, FLEXIBLE AND 
RIGID BEARINGS, COUPLINGS AND GYROSCOPIC 
EFFECTS. CALCULATES CRITICAL SPEEDS AND 
WRITES CONVERGED VALUE. CALLS 
SUBROUTINE GRA AND PPLOT FOR APPROPRIATE 
CURVES AND DRAWINGS. 

£*********************************************************************** 
REAL ID ( 200 ) , MASS(200) 

INTEGER STA(IO) 

EXTERNAL SU8R 
COMMON IO,MASS,STA 
DIMENSION XA ( 500 ) 

COMMON RAT (10),CS(3), DX( 200 ) , OD! 200 ) , SSI 200 ) , E ( 200 ) , RHO! 200 ) , 

IT AOD (200)»TAID( 200 ) » ARE A! 200 ),AI!2G0), DEL 1200), BRG ( 10) ,BC!10), 
2ZZI200) ,EI i 200) ,GAM (200) , LOCI 200 ) , THIC < 200 ) ,DEN( 200) ,0DIA(200) , 

3 DIA( 200) ,CRIT{ 6,3) ,K0 ( 10 ) , DY ( 58 ) , X ( 500 ) , Y1C 500 ) , Y2 ( 500 ) , KL ( 200 ) , 
4EM! 400 ),YP!10),SL! 200 ) , SLL ( 200 ) ,0N !10)*TITLE!12) 

COMMON PI , KEY, NNN, IC, SCALE, ICM, ISK, IJJ,KA, INN, IXY,RPM , NO,NUM, 
1TP,RP,KB, WW, IR, J ,KK,DRPM,FF,L , TOTAL, I XZ, NOO, RPMF , MOD 
COMMON/CLPLOT/XPEN, YPEN»NX,NY, I PEN, XL ABE L( 10) ,YLAB£L( 10) 

400 FORMAT ! 1H0,26HA DISCONTINUITY OCCURS AT ,F10.0,33HIN THE EXCESS MO 
1MENT VS RPM CUR VE/ 1H0 , 27HDATA FOR THIS CURVE FOLLOWS) 

501 FORMAT ( 1H0 , 75X, 20H1ST CRITICAL SPEED -F8.0) 

502 FORMAT ( 1H0, 75X, 20H2ND CRITICAL SPEED * F8.0) 

503 FORMAT ( 1H0 , 75X, 20H3RD CRITICAL SPEED * F8.0) 

504 FORMAT! 1H0 , 15X, 22HN0 CRITICAL SPEED FR0M,F8.0, 2HTO,F8.0) 

506 FORMAT! 1H0,16HDID NOT CONVERGE) 

4000 FORMAT !1X,4!F8.0,E15.8) ) 

100 WW=!2.*RP*PI /60.)**2 
DO 99 1=3,10 
99 DY ! I )=BC ( I ) 

I R=1 

J=1 

KK=1 

DY ! 1 ) =0 • 

X ! J) =0 • 

Yl! J )=DY!6) 

Y2! J ) *DY ! 10 ) 

IM=1 

IU=1 

DO 50 K= 1 , L 
OD! 200 ) =X ( J ) 

KK=K 

DY! 2 ) =DX !K)/DEL!K) 

I T=DEL ( K ) 

DO 90 IVP=1, IT 
90 CALL RUNGEK! 8»DY,$UBR) 

I F! SS ( K ) • EQ.O. ) GO TO 44 
DY!3 )=SS(K)*WW*DY!6 )+DY! 3) 

DY!7 ) =SS!K)*WW*OY! 10)+DY!7) 

44 IFtSTA! IU) .EQ.KK ) GO TO 45 
GO TO 96 

45 IF(BRG! IU))81,82,83 
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83 DY(3)=DY(3) -BRG ( IU ) *DY ( 6 ) 

DY ( 7 ) - DY ( 7 ) -BRG ( IU ) *DY ( 10) 

GO TO 95 

82 DY(9)=DY( 9 l-DYC 10 » *DY ( 5 > /DYC 6 ) 

D Y{ 8 ) =0Y ( 8 ) -DYC 4 ) *DY( 10)/DY< 6) 

RAT ( IR)=DY( 10)/DY(6) 

K0( I R ) ~ J 
I R=l R+l 
Y2C J)=0. 

Y1 ( J ) -0 • 

DY(10)*0. 

DY ( 7 ) =0 • 

D Y( 3 ) = . 0001 
DY ( A ) -0 • 

DY(5)-0. 

DY( 6 ) *0 • 

GO TO 95 

81 DYC 7 )*DY ( 7 ) -DY ( 3 ) *DY( 8 ) /DY ( 4 ) 

DYC 10 ) = DY ( 10 )— DY ( 6 )*DY( 8)/DY(4) 

RAT ( IR)=*DY(8) /D Y C 4) 

K0( I R ) = J 
I R* I R+l 
Y2C J ) *DY CIO) 

Y1 ( J ) *0 • 

DYC 8 ) =0. 

DY ( 9 ) *0* 

0Y(3)*0. 

DY ( 4 ) *0 . 

DY(5)*.0001 
OY ( 6 ) =0 • 

95 IU=IU+1 

96 IF(LOCCIM).EQ.KK) GO TO 97 
GO TO 50 

97 DY(4)-ZZ(K)*DY( 5)*WW+DY(4) 

DYC 8 ) -ZZ ( K ) *DY( 9 ) *WW+DY ( 8 ) 

I M* I M + l 

50 CONTINUE 

RAT( IR)»DY(7)/DY(3) 

K0( I R )* J 

IFCIXY.EQ.l) GO TO 295 

EMC INN)=DY(3)*DY(8)-DY(4)*DY(7) 

IFCINN.LT. 2) GO TO 299 
IF (EMC INN)*EM( INN-1 )) 298,299, 299 
298 SO*RP 

RP* ( EM ( INN ) *DRPM ) / ( EM( INN-1 )-EM( INN) l+RP 
W3=RP 
JKK*0 
I XY» l 

W1»S0— DRPM 
W2*EM( INN-l) 

W5*S0 
W6=EM( INN) 

GO TO 100 

295 W4*DY (3 ) ♦DY ( 8 )-DY ( 4 ) ♦DY ( 7 ) 

IF( ABS( W4).GE.A8S(EM( INN) ) .OR. ABS ( W4 ) .GE. ABS4 EMC INN-1) ) ) GO TO 300 
IFCJKK.EQ. 35) GO TO 2968 
JKK* JKK+1 

I F C W2*W4 ) 2200,2200,2201 
2200 W5=W3 
W6=W4 
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GO TO 2202 

2201 W1=W3 
W2*W4 

2202 PRP* W3 

RP*— W6* C Wi— W5 ) / (W2— W6)+W5 
W3*RP 

I F( ABS { PRP— RP ) • LE *9 • ) GO TO 2967 
GO TO 100 

2967 CRIT(NNN, IJJ)*PRP 
I XY*0 

RP * SO 
GO TO 803 

2968 WRITE (6.506) 

I XY*0 

GO TO 39 

803 CS(IJJ)*CRIT (NNN* I JJ) 

KEY*0 
CALL GRA 
705 IJJ*IJJ+i 

I F C I JJ.EQ.4) GO TO 78 
299 INN*INN+i 
RP*R P+DRPM 
N0*N0-1 

IF(NO.NE.O) GO TO 100 
78 INN*INN+1 
OCN*NOO 

ZN* OON*DRPM +RPM 
IFCCSm.EQ* 0.) GO TO 666 


WRITE<6, 

5011 

csm 



IF(CS(2) 

• EQ. 

0. ) GO 

TO 

667 

WRITE (6 * 

5021 

CS (2 ) 



iFicsm 

• EQ. 

0. ) GO 

TO 

667 

WRITE ( 6 » 

503) 

CS(3) 




667 GO TO 39 

666 WRITE (6*504) RPM.ZN 

39 KEY*1 

CALL PPLOT 
CALL GRA 
NP*INN-1 
I CU* INN-2 

I F( 1 XZ • EQ* 0) GO TO 451 
XA( 1 )*RPM 
DO 25 1*1. ICU 

25 XA(I+1)*XA(I)+DRPM 

WRITE (6 f 4000 ) (XA(I), EM(l), 1*1. NP) 
451 IFINNN.EQ.6) GO TO 900 
RETURN 
900 KEY*— 1 

CALL GRA 
RETURN 

300 WRITE (6,400) SO 
IXZ*1 
I XY*0 
RP*SO 
GO TO 299 
ENO 
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noon 


C*«t4#*4l******4 ***************»***********»4‘** ************* + **4t*** **♦*♦* 

SUBROUTINE GRA 

PREPARES DATA FOR THE DEFLECTION 
ANC BEARING STIFFNESS VS CRITICAL 
SPEEDS CLRVES 

Q* 4* * 44 * 4 *** 4 * **4***************4* 4**4 **************** ****************** 

SUBROUTINE GRA 

REAL ID ( 2GC ) * MASSI200) 

INTEGER STA(IO) 

COPPCN IC *PASS*STA 

COPPCN RAT ( 10 ) ,CS ( 3 ) ♦ DX (200 ) *00 (200 ) * SS ( 200 ) * E ( 200 ) » RHOt 200 ) , 
1TA0D I 200 ) * TAI D (200 ) y ARE A (200 ) * AI (200) * DELI 2001 *BRG< 10),BC( 10), 

2 ZZ < 200 ) , E I < 200 ) , G A M < 200 ) * LOC ( 200 ) * TH I C C 200 ) » DEN ( 200 > * OD I A ( 200 ) t 

3 0IAI200) ,CRIT(6,3) ,K01 10) * DYC 58 ) , X 1 500) *Y1 1 500 ) * Y2I 500 ) *KL < 200 ) , 
4EIM4C0) ,\P(10) * SL < 200 ) * SLL < 200) » ON ( 10 ) * T ITLEl 12 ) 

COPPON PI *KEY*NNN*IC, SCALE* ICM, ISK, I JJ,KA, INN, IXY,RPP » NO*NUM, 

1TP *RP ,KB * WW,IR* J,KK*DRPM,FF,L * TOT AL , I X Z , NOO , RPMF * MOO 
COPPCN/SPEC L /TEST, ORGSET, SPA SET 

COPPCN/CLPLOT/XPEN*YPEN,NX,NY,IPEN,XLABEL< 10) ,YLABEL( 10) 

01 PENS ION BELX(IO) , BELY (10) *CA1 18 ) *KKK(40 ),P< 20) ,DEF< 500 > 

1 * A Z ( 500 ) *ABELX<10) , ABEL Y (10) 


DATA ( 

ABELX1K), K*1 ,10) /6H 


*6H 


* 6H 

,6H 

* 

1 6HLENGTF.6H. *6H 

* 6H 


* 6H 


,6H 

/ 

DATA I 

A6£LY(L)*L*1, 10) / 6H 


*6H 


*6H 

* 6H 

* 

1 6HDEFLEC *6HTI0N* *6H 

* 6H 


* 6H 


* 6H 

/ 

DATA 1 

eELY(L)*L*l,10)/6H 


*6H 


* 6H 

,6HCRITIC, 

1 6HAL 

* 6HSPEEDS *6H 

*6H 


* 6H 


* 6H 

/ 

DATA ( 

BELX ( K ) *K*1 * 10) / 6H 


,6H 


* 6H 

* 6HBEAR IN* 

1 6HG 

*6HSTIFFN,6HES$ 

*6H 


, 6H 


*6H 

/ 


IF (KEY) 55 *90*777 

59 IF (CRITIItl).EQ. 0.)G0 TO 777 
DO 100 1*1*6 
AZ(I)*ON(I) 

AZtI+6)*AZ( I ) 

IOC AZ(I412)*AZU) 

11*1 

DO 76 1*1*10 
XLABELU )*BELX(I ) 

76 YLABEL(I)*BELY(I) 

DO 20 KV* 1*3 

DO 2C 1*1*6 

OA < 1 1 ) *C R I T ( I * K V ) 

IF < C A ( 1 1 ) • N E • 0 • ) GO TO 20 
CA< II )*CA< II-l) 

AZ C 1 1 ) *A2 ( 1 1-1 ) 

20 11*11+1 
KKKC1 )=1 
KKK < 2 ) *1 
KKK < 3 ) =3 
KKK < 4 ) *1 
KKK ( 5 ) =1 
KKK (6 )*6 
KKK 17) *6 
KKK C 8 ) *6 
P < 1 ) *3« 

PI2)*10. 

PI 5 )*10* 

p< e)*io. 

DO 77 1*5*13 
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7? P(I)*0. 

PI 141*90. 

XPEN=0. 

YPEN*0. 

NX*-60 
NY*8C 
FF*2 • 

GRGSET=0 . 

I F { I CM* EC • 1) GO TO 777 
CALL CALFLTIAZ,CA,KKK,P) 

GO TO 777 

90 DO 91 1=1, NOP 

91 AZ(I )*XC I)*SCALE 
NP*1 

LZ*1 

N^KC 1 LZ 1 

79 DO 88 I=KP,N 
88 DEF( n=Y2(I)-Yim*RAT(LZ) 
LZ=LZ+1 

IF IKC I LZ 1 .EQ.O) GO TO 74 

NP*N*1 

NsfKC(LZ) 

GO TO 79 
74 TE $ T = 50. 

0RGSET=1 
DO fiCO 1*1,10 
XLA8ELI I 1*ABELX( I ) 

60C YLABEL II1=A8ELYCI1 
KKK ( 1 1=4 
KKK (21=1 
KKK (31=1 
KKK C 4 1 =1 
KKK I 5 1 = I J J 
KKK ( 6 ) = J 
Pf 11=1. 

PI2 )=TP*SCALE 
PI3)*0. 

PC41*TP*5CALE 

P 1 5 ) * 10. 

PI6>*-5. 

PI 7 1 = 5. 

PI8M10. 

DO 877 1*9,13 
877 Pfl>*0. 

PI 14 1 =90. 

XPEN*0. 

YPEN*0. 

NX*-60 
NY*6C 
H I =C • 

FF = 1 . 

DO 333 1*1 , J 

333 HI*AKAX1 (AESIDEFII 1 1 ,HI 1 
DC 334 1*1, J 

334 DEFI I1*5./HI *DEF 1 1 1 

IF I ICM.EC. 1) GO TO 777 
CALL CALFLT IAZ,DEF*KKK*P) 
777 RETLRN 
END 



nooon oooooooo 


c ***,**i:**S*Iiil*S*i2I****************************************»********* 

SUBROUTINE PPLQT 
PREPARES DATA FOR SHAFT DRAWING 

CALPLT IS A GENERAL PURPOSE ROUTINE FOR PRODUCING PLOTS ON THE 
CALCCMP MODEL 670 DIGITAL PLOTTER* 

CALL C ALF LC T ( XDOWN * Y ACROS » KKK ♦ P ) AR(fKa 

WHERE — XDOWN IS AN ARRAY CONTAINING THE DATA FOR THE ABSCISSA 
YACROS IS AN ARRAY CONTAINING THE DATA FOR THE ORDINATE 
KKK AND P ARE ARRAYS THAT CONTAIN PLOTTING INFORMATION 

THERE ARE TWO COMMON BLOCKS 
COMMON /CL PLOT / 

CCMMCN/SPECL/ 

PLC T * SYMBOL t NUMBER t AND AXIS ARE SUBROUTINES USED BY CALPLT 
4*4 4 #***********: 


c*** 


INTEGER STA(IO) 

REAL I D ( 2CC ) ♦ MA$S(200) 

COMMON RAT ( 10) *CS(3) « DX (200 > 1 00 (200 ) » SS ( 200 ) t 6 ( 200 ) t RHO( 200 ) t 
1TA0C ( 200 )fTAIO(200)*AREA(200)tAI (200 ) ♦ DEL < 200 ) • BRGC 10 I » BC( 10 )♦ 
2ZZ(2CO)»£I(2QO) » GAM (200 ) *LOC (200 ) * THIC ( 200 ) » DEN ( 200 ) »ODI A( 200) * 

3 DIA(200t,CRIT(6,3),K0ll0),DV(58),X(500>.Vl(500),V2(500),KL(200), 

4EM(4CO)*VP(lO)tSL(2O0)*SLL(2OO)»ON(lO)*TITLE(12) 

COMMON P I f KEY tNNNf 1C » SCALE * I CMt I SK 1 1 J J »KA » INNf I XYtRPM » NO»NUM* 
lTP*RPtKBt WW» IR * J»KK*DRPM f FFtL t TOTAL t IX Z *N00 » RPMF t MOD 
COMMON/SPECL /TEST *ORGSET *SPASET 

CCMMCN/CLPLOT/XPEN»YPENtNX»NY* I PEN f XL ABEL ( 10) *YLABEL ( 10) 

DIMENSION 00(400 ) »ODL (1600) * Z(200 ) » Z0( 400) *KKK( 40 ) ,P ( 20) 
IF(ICM.NE.C) GO TO 38 


FF*C • 

TEST*50. 

CRGSET=1* 


JJJP*1 

DC 120 I * 1 *KB 
JC*ST A ( I ) 

TCT=C. 

DC 131 KF P *1 * JC 
131 TGT*CX(KFP)*SCALE+TGT 
YP ( J J JP ) * TOT 
12C JJJP=JJJP+1 
Jf>*l 

DO 5C1 I J-l *L 
CD < JP ) = ID ( I J) 

0DL( JP)®CD(IJ)**5*SCALE 
JP*JP+1 

IF (TAID( I J).EQ.O*) GC TO 80 
CD( JP) =TA ID ( I J) 

GO TC 40 

80 CO ( JP ) - It ( I J ) 

AO IF (TAOD( I J ) .EQ.O. ) GO TO 81 
00L( JP)*1AC0(I J)**5*SCALE 
GO TC 501 

81 ODL ( JP )*CDL ( JP-1 ) 

501 JP*JP+1 

DO 88 l A* 1 *L 
68 Z( IA )=0X( I A )*SCALE 

ZC(1)*0. 


JP *2 

DO 48 I A* 1 » L 

ZQ ( JP ) =ZC ( JP-1)+Z(IA) 
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noon 


48 


JP=JP+1 

ZC< JP)=ZC< JP-l) 

JP*JP+1 
KH=2*L 

OG 66 I A * 1 , KH 
KM*KH+I A 
66 00 L ( KM ) =C 0 ( ! A ) *. 5*SC ALE 

KMM=KM+i 
DG IQ 8 I 4* 1 *KH 
ODUKMM)*-CDLtIA> 

106 KMM*KMM+1 
KWK ( 1 ) =4 
KKK ( 2 ) =0 
KKK ( 3 ) =4 
KKK I 4 ) = 1 
KKK(6)*KF 
Pll)=l. 

P«2)»TP*SCALE 
P ( 3 ) =0. 

P C4 ) =TP*SC ALE 

Pt5)=10.C 

P(6)=-5. 

PI7M5. 

PC6 )=10.C 
DO 771=9,13 
77 P<I)=0. 

Pi 14)=90. 

XPEN=0*G 
YPEN*0.0 
NX*-60 
NY = 6C 

CALL C ALP LT ( ZO ,CDL t KKK* P ) 

CALL CALTIT 
RG*1. /SCALE 
P< 6)^-5. /SCALE 
P|7>=5. /SCALE 

CALL AXIS tP<9) ,P(10> ,XLABEL,NX,P(2>, P< 13 > , P ( 3 ) , RG, P < 8 ) ) 
RG*(P(7)-P(6> >/P<5> 

CALL AXIS i Pill) tPi 12) ,YLABEL,NY,P(5) , P { 14 ) * P ( 6 ) , RG, P < 8 ) ) 
CALL P L0T(P<2)+10.,0.,-3> 

38 RETURN 
END 


c********* ***************** ************************************ ********* 

SUBROUTINE SUBR 

THIS SUBROUTINE IS CALLED BY RUNGEK 
TO CALCULATE VALUES OF THE DERIVATIVES 
FOR THE EQUATIONS OF BEAM FLEXURE. 
£*********************************************************************** 
REAL I D( 200 ) , MASS(200) 

INTEGER STAUOI 
COMMON ID, MASS, STA 

COMMON RAT ( 10 ) , CS ( 3 ) , DX( 2001 , OD C 200 1 , SSI 200 >, E { 200 > ,RHO ( 200 ) , 
1TA0D(200) ,TAIDC 200 >, AREA! 200), AI< 200 ),DEL( 200), BRG( 101 ,BC(10), 
2ZZ(200),EI( 200 ) , GAM ( 200 ) , LOCI 200 > , THIC ( 200 ) ,DEN { 200 ) ,ODI A ( 200) * 
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80 


3 DI A ( 200 ) * CRIT ( 
4E M ( 4-00 ) * YP ( 10 
COMMON PI, KEY 
1TP,RP * KB , WW 



KK~KK 

IF(KA.EQ. I) GO TO 80 

AA=OD(KK )+SL(KK ) *( DY( 1 ) -00(200 ) ) 

EE* I D(KK)+SLL(KK )*( DY ( 1 )-0D( 200) ) 

G*3» *E ( KK ) /8. 

TZ* ( EE/ AA ) **2 

A2*AA*AA 

D2*EE*EE 

CROSS=( A2-D2)*PI/4. 

F= ( 9 . * { 1 . +TZ ) **2+24.*TZ )/(8**(l.+TZ)**2) 
GAM(KK)*F/ (G*CROSS) 


A4=A2*A2 


D4*D2*D2 

MASS(KK )*( A2-D2 ) *P I *RHO (KK ) /4. 
El ( KK ) * ( A4-D4 )*PI/64.*E (KK ) 
DY(11)*MASS(KK)*DY(6)*KW 

DY (12 )*DY ( 3 ) 

DY(13)*DY(4)/El (KK) 

DY( 14 ) *DY ( 5 ) — DY ( 3 ) *GAM ( KK ) 

DY( 1 5 ) —MASS ( KK) 4DY ( 10 ) *WW 
DY (16 ) *DY ( 7 ) 

0Y(17)*DY(8)/EI (KK) 

DY(18 )*DY(9)-DY( 7)*GAM(KK> 

RETURN 

END 


* 


SUBROUTINE RUNGEK (NN , Y , SU0R ) ************ *** 

C************* 4 ************************* ******************************** 

SUBROUTINE RUNGEK 

PERFORMS THE NUMERICAL INTEGRATION 
THE DIMENSION OF Y IN MAIN MUST BE AT LEAST 7*N+2 
Y(1)*X, Y ( 2 ) =DX , Y(3) , . • . , Y ( N+2 )*Y1 , • * • , YN 

C CAtL^UBlTlI^MM SET UP THE INITIAL FIRST DERIVATIVES 

^♦♦IjJiU****************'************************************************ 

REAL ID(2CC), MASS(200) 

INTEGER 5 T A ( 10 ) 

COMMON RAT(IO) ,CS(3) , DX ( 200 ) ,0D ( 200 ) , SS ( 200 ) , E ( 200 ) , RHO ( 200 ) t 
IT AGO ( 200 ) ,TAID(200) , AREA (200) , A I ( 200 ) , DEL ( 200 ),BRG(10),8C(10), 

2ZZ (2C0) ,E I (200) ,GAM(200) ,LOC ( 200 ) , TH IC ( 200 ) , DEN < 200 ) , ODI A ( 200 ) , 

3 0 I A (200 ) »CRIT (6*3) ,KC(10) , D Y ( 58 ) , X ( 500 ) ♦ Y 1 ( 500 ) , Y2 ( 500 ) ,KL ( 200 ) , 
afm (400) ,YP(10) ,SL(20C) ,SLL(200) , ON (10), TITLE! 12) 

COMMON PI , KEY, NNN,IC, SCALE, ICM, I SK , I JJ,KA, INN, IXY,RPM , NO,NUM, 
1TP*RP,KB, WW , I R , J , KK,DRPM, FF , L , TOT AL , IXZ , NOO , RPMF , MOD 

DIMENSION Y ( 1 ) 

CALL SUSP 


N=?NN 

DC 1C I *1 # N 
I l=6*N+2* I 
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1C Y«Il)*YU + 2) 

00 12 1*1 tN 

I1*N42+I 
I2*2*N+2+I 
12 YU2)-Y(2)*Y(Il) 

Yll)*Y(l)+Y(2)/2. 

DC 1 A 1*1, N 

tl*2*N+24l 

E2*6*N+2*I 

1 A YU + 2)=YU2) + .5*YU1) 

CALL SUBP 
00 16 1*1 1 N 
Il*N*24l 
I2*3*N+24I 
16 YII2)*Y<2)*YCI1) 

00 18 1*1, N 

Il=3AN424l 

I2*6*N424l 

18 YU42)*YII2) + .5*YCU) 

CALL SUBP 
00 2C 1*1, N 
I1*N42+I 
I2*A*N+24 I 
2C Y« 12 )*Y(2 )*YC 11) 

YC1 )*Y<l) + Y(2)/2* 

DO 22 I = 1 , N 

ll*4*N+2+l 

I2*6*N+24I 

22 Y 1 1 + 2 1 *Y (I2)4Y(I1) 

CALL SUBP 
DC 2A 1 = 1 ,N 
I1*N424I 
I 2*5*N+24 I 
2 A YC 12 ) *Y ( 2 ) * Y ( I 1 ) 

DO 26 I = 1 , K 

Ii*2AN+24l 

I2*3*N424l 

I3*A*N4-24 I 

IA»5*N424I 

I5*6*N+24l 

26 Y I l42)*YU5) + m I1)+2.*Y(I2)42.*Y( I3)+Y( IA) 1/6. 
J*J4l 

Yl< J)=Y(6) 

Y2(J)*YU0) 

XI J)*Y(11 

RETURN 

END 



SUBROUTINE CALTIT 

C **4 4**4|*****44 ♦****♦**♦*********♦♦♦♦♦ *****♦*♦♦♦♦*********♦♦**♦**♦*♦♦♦** 
C SUBROUTINE CALTIT 

C WRITES HCLLERITH INFORMATION ANO NUMBERS 

C ON PLOTS* 

4 4****4«* *4 4 ♦♦♦4t*****»***4t**4<*******4t**4«4t**4>*4t**4**4»**4t *♦***♦♦*♦* ***4t 

REAL I D ( 2C0 ) y MASS(20O) 

INTEGER !TA(10) 

COMMON ICtMASSfSTA 

COMMON R*T (10) ,CS<3> , OX 1200 ) 9 00 1 200 ) 9 SS ( 200 ) 9 E ( 200 ) 9 RHO ( 200 ) 9 
1 TACO < 200 ) 9 TAI D (200 ) 9 AREA (200 ) 9 A I (200 ) 9 0EL( 200 ) 9 BRG( 10 ) 9BC( 10 ) 9 
222 ( 2G0) ,E I (200 ) 9 GAM ( 200) 9 LOC (200 ) 9THIC( 200 ) ,OEN( 200 ) t ODI A( 200 ) 9 
3 0 1 A ( 200 )9CRIT(6 f 3)9K0(10)9DY(58) ,X( 500 I 9 YII 500 > ,Y2( 500) 9 KL ( 200 ) 9 
A EM ( A 00 ) t>P (10) ,SL<200 ) f SLL(200) 9 ON (10) f T 1TLE ( 12 ) 

COMMCN PI ,KEY f NNN 9 IC 9 SCALE 9 ICM t ISK 9 IJJ 9 KA 9 INN f IXY 9 RPM 9 NO 9 NUM 9 
ITP 9 RP 9 KB 9 WW 9 IR 9 J 9 KK,DRPM,FF t L , TOT AL , IX Z , NOO , RPMF , MOD 
IF ( FF. EQ. 2. ) GO TO 155 
IF ( FF.EQ. 1* ) GO TO 177 
176 CALL SYMBOL ( 1. ,9.5,. 1 ,TI TLE,0. ,72) 

XXX*C. 

OG 50 I * 1 9 L 

X*X*CX( I >4<SCALE*XXX 

IF ( SS ( I ) .EC.O. ) GO TO 50 

Q^SSII) 

CALL SYMBOL ( XXX,0. , • 1 9 1 1 HEXTRA M ASS* # 90. 9 1 1 ) 

CALL NUMBER (XXX 9 i. 29 .lf Q t90. 96 ) 

50 CONTINUE 
YY*9.3 

EF(CSd) .EC.O. ) GO TC 40 
KVV*I JJ-1 
00 1C I *1 9 KV V 

CALL SYMBOL (1. ,YY *.1,1 5H CRITICAL SPEED* f 0. 9 16 ) 

II*I 

£Q*CS < I ) 

CALL SYMB 0 L( 3 . 3 # YY 9 . 1 rII f0.9“l) 

CALL NUMBER(3.59YY,.l9QQtO. 9 O) 

10 YY*YY— .3 

AO CALL SYME0L(1. 9 YY, . 1 9 11HT0TAL MASS*90.9ll) 

CALL NUMBER(3.5,YY,.l4T0TAL90. 96 ) 

IF(ISK .EQ. 0) GO TC 35 
YY*YY-.3 

CALL SYMeOL( 1 . 9 YY 9 . 1 f 49 HGYROSCOPlC EFFECT— SEE PRINTED OUTPUT FOR 
1DET A I LS 9 C . ,49 ) 

35 IF(KB.EQ.O) GO TO 177 
DO 17 1*1, KB 
IF (BRG ( I ) ) 16919,18 

18 YPP*YP ( I ) 

CALL SYMBOL (YPP,5. , . 1 ,8HBEAR ING* ,90. ,9 I 
Q 1 C*BRG( I ) 

CALL NUMBER ( YPP , 6 . 0 , . 1 9 QQ ,90. 9 0 ) 

GO TC 17 

19 YPP*YP ( I ) ^ .... 

CALL SYMBOL (YPP,5* ,.l ,13HRIG ID BEARING,90. 9 13) 

GG TC 17 

16 YFP*YP(I) 

CALL SYMBOL (YPP, 5. 1 ,8HC0UPLING ,90. * 8 ) 

17 CONTINUE 
GO TC 7 

177 CALL SYMEOLd. , 9. 5 , . 1 ,T I TLE ,0. ,72 ) 
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7 RETURN 

155 CALL SYMBOL <1. ,9.5 ,.l ,TITLE,0. ,72) 

CALL SYMBOL (l.»9.2».l,l,0.,— 1) 

CALL SYMBOL (1. 2,9.2, .1 , 20HF I RST CRITICAL SPEED, 0., 21) 
CALL SYMBOL U., 8. 9,. 1,2,0. ,-11 

CALL SYMB0L(1.2,8.9,. 1 ,2 1HSEC0ND CRITICAL SPEED, 0., 22) 
CALL SYMB0U1.,8.6,.1,3,0.,-1) 

CALL SYMBOL 11. 2,8.6, .1,20HTHIRD CRITICAL SPEED, 0., 21) 

GO TO 7 

END 
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